function sandbox = fsb_calculate_full_map(sandbox,idat,hrf_pred,statsel)

% FSB: calculate a statistical map for the whole volume
%
% EXAMPLE:
% sandbox = fsb_calculate_full_map(sandbox,idat,hrf_pred,statsel)
%
% INPUT
% sandbox: sandbox experiment struct
% idat: 4D-Image data
% hrf_pred: hemodynamic predictor number selected
% statsel: method for map calculation
%
% OUTPUT:
% sandbox: sandbox struct containing calculated map
%
% CALLED BY:
% FSB.m
%
% NOTES
% Use this function to calculate a full map instead of only visible sections
% it takes a little longer to compute but speeds up subsequent processing
% Three options are available for computation which are quite different
% in processing time. Correlation is fastest (statsel == corrc), the two
% other options (statsel == robreg || statsel == glm) are quite slow.
%
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
%
% $ Revision 1.0
%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

h = waitbar(0,'Calculating statistics for whole volume...');
tic ;
sandbox.hemodynamics = single(sandbox.hemodynamics);
%idat_single = single(idat);
[a b c d] = size(idat);

switch statsel

    case 'corrc';
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        % Correlation coefficient (fastest)
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

        if d>5 % Use this code path for larger datasets

            waitbar(1/10)
            idat2 = shiftdim(idat,3);
            waitbar(2/10)
            idat3 = reshape(idat2,d,[]);
            waitbar(3/10)
            test = fsb_corr2d(sandbox.hemodynamics,single(idat3));     
            %test = corr(sandbox.hemodynamics,idat3,'rows','complete');
            waitbar(7/10)
            test2 = reshape (test,[],a,b,c);
            waitbar(8/10)
            test3 = shiftdim(test2,1);
            waitbar(9/10)
            test3= squeeze(test3);
            sandbox.stats.r_cor =test3;

        else % Use this code path for smaller datasets

            for hem = 1:size(sandbox.hemodynamics,2)
                waitbar(hem/size(sandbox.hemodynamics,2))

                for y=1:b;
                    for x=1:a;
                        for z =1:c;
                            if sandbox.stats.r_cor(x,y,z,hrf_pred)==0;
                                R = corrcoef(squeeze(single(idat(x,y,z,:))),sandbox.hemodynamics(:,hem));
                                if abs(R(2,1)) <0.01;
                                    R(2,1) = NaN;
                                end
                                sandbox.stats.r_cor(x,y,z,hem) = R(2,1);
                            end
                        end
                    end
                end
            end
        end

        close(h);


    case 'glm';
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        % GLM: slow
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        for z =1:c;
            waitbar(z/c)
            for y=1:b;
                for x=1:a;
                    if sandbox.stats.glm_prob(x,y,z,hrf_pred)==0;
                        [besti2,dev,statsi] = glmfit(sandbox.hemodynamics,squeeze(single(idat(x,y,z,:))));
                        sandbox.stats.glm_prob(x,y,z,:) = statsi.p(2:end);
                        sandbox.stats.b_est(x,y,z,:) = besti2(2:end);
                    end
                end
            end
        end

        close(h)


    case 'robreg';
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        % robust regression: slow
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        for z =1:c;
            waitbar(z/c)
            for y=1:b;
                for x=1:a;
                    if sandbox.stats.rob_prob(x,y,z,hrf_pred)==0;
                        [robr,statsi] = robustfit(sandbox.hemodynamics,squeeze(single(idat(x,y,z,:))));

                        if statsi.p(2:end) ==0;
                            statsi.p(2:end) = NaN;
                        end

                        sandbox.stats.rob_est(x,y,z,:)=robr(2:end);
                        sandbox.stats.rob_prob(x,y,z,:) = statsi.p(2:end);
                    end
                end
            end
        end

        close(h)

end
t = toc;
disp(['Full map calculation done in ' num2str(t) ' seconds']);

end
